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Abstract 

Recently, it has been proven [R. Soc. Open Sci. 1 (2014) 140124] that 
the continuous wavelet transform with non-admissible kernels (approximate 
wavelets) allows for an existence of the exact inverse transform. Here we 
consider the computational possibility for the realization of this approach. 
We provide modified simpler explanation of the reconstruction formula, re¬ 
stricted on the practical case of real valued finite (or periodic/periodized) 
samples and the standard (restricted) Morlet wavelet as a practically im¬ 
portant example of an approximate wavelet. The provided examples of ap¬ 
plications includes the test function and the non-stationary electro-physical 
signals arising in the problem of neuroscience. 
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1. Introduction 


The continuous wavelet transform (CWT) with the standard (restricted) 
Morlet wavelet 0(£) — exp — £ 2 /2) is defined as the integral 
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where a and b are the scale and the shift correspondingly, ojq is the central 
frequency. 

It is one of the most powerful modern tools of signal processing espe¬ 
cially adjusted to the extraction of instant oscillating patterns [lj, 2fl since 
the transform (JT]) of the harmonic oscillation f(t ) = exp (hut) results in the 
complex function 

, (au— tun)^ 

w(a, b) = e lub e~ 2 

The modulus of this function has a maximum, which allows for the determin¬ 
ing of signal’s frequency lu = u} 0 /a max and the phase coincides with signal’s 
one. 

The modern applications of the continuous wavelet transform are focused, 
in particular, on a study of environmental time series [§, dj, geo- and astro¬ 
physics [5|,l6i,[7|, biophysics [8j,|9, 10] and neuroscience, see an extensive review 
in the recently published book jll| . 

At the same time, the actual problem is not restricted by the search of 
oscillating patterns localizations. It is important to extract revealed struc¬ 
tures from the b ackg round consisting a noise, global oscillations and inho¬ 
mogeneities, etc. 12, 13, 14, lil . 

However, the conventional approaches to the inversion of the wavelet 
transform with the standard Morlet wavelet have some principal difficulties 
from the point of view of functional analysis. To be applicable in a classi¬ 
cal inversion formula, a wavelet function -0 should satisfy the admissibility 
condition [lj 
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where the asterisk denotes a complex conjugation. Then the classical inver¬ 
sion of the wavelet transform W^f = w(a,b) with a wavelet function ifj is 
written as 



However, the integral diverges for the standard (restricted) Morlet wavelet. 
On the other hand, the alternative inversion formula for the CWT 




is proven in [16] under a mild natural conditions on a wavelet function ip. 

The principal aim of this paper is to show how this approach can be 
realized in applications. For this reason, we adapt the proof of (J2J) to the 
case of real-valued functions with a finite support in such a way that its 
line of reasoning and the result can be straightforwardly used as a practical 
algorithm for the reconstruction of a function from its wavelet transform with 
the Morlet wavelet. 

The paper is organized as follows. Section 2 shortly discusses the im¬ 
plementation of the CWT with the Morlet wavelet in the application to 
real-valued functions determined within a finite interval (or periodic on M) 
and then presents the simple proof the exact reconstruction formula that does 
not require the wavelet admissibility condition. Its formulation allows for the 
simple computational realization, which is presented in Section 3. The exam¬ 
ples comprise the test one and the solution of the modern practical problem 
in the held of neuroscience. Neuronal electric activity is characterized by the 
high non-linearity and complexity, extraction of main oscillatory components 
and localization of them is very important for reveal of interaction dynamics 
between different cells in neuronal network Q. The final section gives some 
outlooks for the perspectives for the applications of the proposed method in 
computational physics. 
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2. Direct and inverse continuous wavelet transform of periodic 
functions 

2.1. CWT expansion 

Let us consider a real-valued function f(t) with the zero mean given by 
its expansion into the Fourier series 

OO 

m = E A n cos(co n t) + B n sin (u} n t), (3) 

n=1 


where oj n = ‘Inn/T. 

This function is periodic with the period T. As well, it can be consid¬ 
ered as a function determined within the interval t G [0, T] and periodi¬ 
cally extended over all M. The last point of view is applicative in practical 
computations since one can deal with finite samples only. The assumption 
of the zero mean is also not restrictive for the problems of computational 
physics since it can be easily achieved by the exclusion of the averaged value: 

m - t-' j 0 T f(t)dt -> /(f). 

For the wavelet analysis of spectral components, it is most convenient to 
operate with the complex analytical counterpart of the function f(t) obtained 
via the Hilbert transform f a (t) = f(t) + iH [/(£)]. 

Since H [cos(x)] = sin(x), if [sin(x)] = — cos(x), the Hilbert transform 
applied to Eq. ()3]), accompanied with Euler’s formula, gives 

OO OO 

fait) = Y J {A n - iB n ) e iUnt = J2 C n e^\ (4) 

71=1 71=1 

Note also that Euler’s formula applied directly to Eq. d3[) results in the 
representation 

OO 

m = r E C n^"‘ + CJe-”" 1 . 

71=1 

where asterisks denotes the complex conjugation. 

Thus, one does not need in practice to apply any special additional proce¬ 
dure for the Hilbert transform of a given sample. It is enough to evaluate the 
discrete Fourier transform (say, via FFT algorithm) and to cut off elements 
corresponding to negative frequencies. 
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The substitution of the series flD) into the integral ([[]), 
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after close form computation of Poisson’s integral f+°° exp(— x 2 )dx = yfn, 
results in the desired formula for the continuous wavelet transform with the 
Morlet wavelet 


w(a,b) = y, C„ 


(cun.a—uj, 
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(5) 


n= 1 


The expression ([5]) is known and usually applied in the numerical algo¬ 
rithms of CWT computations based on the intermediate FFT, e.g. realised 
in WaveLab package [l8j]. However, the preprocessing cut-off procedure, de¬ 
scribed above, allows for the shortening of the time/memory consumption 
because of the twice shortened sample’s length in comparison with the stan¬ 
dard application of Eq. (J5]) to the initial function f(t) directly. At the same 
time, this procedure does not loose any information about f(t) since the 
complex coefficient C n = A n — iB n contains information on both real-valued 
ones. 

As well, the representation ([5]) provides the most clear way to obtain the 
exact inverse continuous wavelet transform. 


2.2. CWT reconstruction 

Now, we aimed to obtain the practically applicable CTW reconstruction 
formula basing on the series (J5]). Its partial differentiation with respect to 
the shift b gives the series 


dw(a , b) 
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Since the included frequencies are strictly positive c u n > 0, the subsequent 
integral over the non-negative scale half-line leads (using Poisson’ integral 
again) to 
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Therefore, taking into account C n = A n — iB n and applying Enler’s for¬ 
mula we get the desired formula for the exact inverse CWT with the Morlet 
wavelet, which does not depend on the admissibility condition: 


fit) 



dw(a, b ) 
db 


da 


( 6 ) 


Note that it reduces computational complexity in comparison with the full 
abstract representation ([2]) since (0 does not require the additional inverse 
Hilbert transform. 

Moreover, Eq. flSJ) allows for the reconstruction of function’s instant fea¬ 
tures revealed by an exploration of the result of the direct CWT. Since the 
partial derivative is a local operator, it is possible to consider w(a,b) with 
b G [bminbmax], the interval, which contains the feature, which we are inter¬ 
ested in. Thus, restricting (jSJ) on the mentioned b, we reconstruct /(f) within 
b G [ tmm = b min tmax = b max ] without boundary disturbances originated from 
the conventional inversion methods based on convolutions. 

As well, since the Gaussian terms in fl5]) are fast decaying functions with 
the maximum in a n = u} 0 /u) n , the restriction of the integral in (EJ) on the in¬ 
terval [a min a max } reconstruct the oscillating components with the frequencies 
from \ujq/ a max loq/ amin\- 

Combining this restriction with the choice of some interval of the shift 
or even considering some non-rectangular region a = a(b), it is possibly to 
reconstruct and analyse very specific details of a studied function /(f). 


3. Tests and applications 

3.1. Test non-stationary oscillations 

As a first numerical example, let us consider the function 


/(f) = e 4t cos(207rf) + yji 2 j (f) sin(407rf), 


(7) 


sampled within the interval f G [0, 1]. 
Here 


X[tb,t e ] {t) 


1, f G [fj, f e ]; 

0, t £ [t b , t e j. 


is the indicator function. Thus, the function ([7]) consists of the decaying oscil¬ 
lations over all studied interval and the stable oscillations of unit amplitude 
over its middle third only. 
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Figure 1: The function 0 and the modulus of its CWT (dark regions correspond to the 
larger values of magnitude). The dashed rectangle marks the maximum line corresponding 
to the localized periodic components with the constant amplitude. 


The function (ED is sampled in 512 = 2 9 equispaced points that allows 
for the direct application of FFT algorithm in further computations, ft 
should be pointed out that the local boundary values of the function (0 are 
sufficiently different that can result in the drastic boundary effects due to a 
periodization. Therefore, we used the algorithm, which expand the interval 
(0, 1] into [—0.5, 1.5] and function’s continuation with boundary reflections 
|2]. The MATLAB code realizing such CWT is presented in Appendix. The 
central frequency is chosen as ujq = 27T, i.e. maxima lines should be located 
along a = 0.05 and a = 0.1. 

Then, the reconstruction formula 0 is applied to this two-dimensional 
wavelet function w(a,b). The partial derivative with respect to the shift 
b is realized as the three-point central difference scheme (except the left 
and the right boundary points, where the forward and the backward two- 
point scheme is exploited). The numerical integration is evaluated via the 
trapezoidal rule. The result of reconstruction for the full shift-scale range is 
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Figure 2: The inverse continuous wavelet transform via the formula ([6]). Upper panel: 
comparison of the reconstruction (solid line) and the original function 0 (dashed line). 
Lower panel: partial reconstructions of oscillating features of the signal 0. 


presented in Fig. [2j One can see an accurate coincidence of the initial and 
the reconstructed signals. 

As the next step, we process the wavelet transformed function w(a,b) 
within intervals that bound the localized oscillating features. This should re¬ 
sult in the decomposition of the function (J7J) into separate instant oscillatory 
components. 

As the first step, we extract dw(a , b)/db for the interval b G [0.3209, 0.6829] 
and integrate this part of the obtained function within the limits a G [0.0240, 0.052], 
This region is bounded by the dashed rectangle in Fig. [0 The result is drawn 
by solid line in Fig. flower panel). Comparing this localized oscillation with 
the almost everywhere constant amplitude with the corresponding compo¬ 
nent of the function (|7j), one can conclude that the proposed methods reaches 
its goal. The pattern detected as a localized spot in Fig. flower panel) can 
be reconstructed with a sufficient accuracy using the simple differentiation 
and integration within the its bounding box. Certainly, it is impossible to 
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avoid boundary disturbances completely, but they are well localized and do 
not corrupt majority of the reconstructed oscillation. 

The rest component is reconstructed by the subtraction of this local com¬ 
ponent from the result of full inverse continuous wavelet transform. It is 
shown as the dashed line in Fig. flower panel). One can see that, except the 
short bursts originated from transient regions, the exponentially decaying 
oscillation is reconstructed as it should be originally. 


3.2. Neuronal oscillations 

As the application of this method, it would be useful to consider an 
example of realistic biological oscillations, in particular, the neuronal network 
dynamics. 

We analyse the complex rhythm, which arises prior to epileptic events, 
so called very fast oscillations (VFO) jl9|. These oscillations appear sponta¬ 
neously and can be revealed after the removal of slow baseline fluctuations. 
Thus, there exists a problem of the reconstruction of various oscillation fea¬ 
tures from the complex signal representing experimental recordings. More¬ 
over, the considered method should allows us to determine an approximate 
time localization for the switching-on and switching-off of oscillations with 
different frequencies and their synchronous states. 

The experimental data are taken from the work [l9[. The potential oscil¬ 
lations are measured in rat’s neocortex: i) the extracellular recordings from 
the neocortex layer (see Fig. 0 upper panel) and ii) the intracellular record¬ 
ings from the pyramidal cell, so-called intrinsically bursting cell, as shown 
in Fig. [U upper panel. All numerical values are presented in dimensionless 
units, which can be rescaled to the experimental values of the time and the 
potential (Fig.l in 19[) as follows: 1 time unit in Figs. [41151 corresponds to 
4.5 ms; 1 voltage unit corresponds to 0.25 mV in Fig. [4] and to 0.5 /rV in 
Fig. H 

The visual exploration of the wavelet modulus plots, Figs 141151 (middle 
panels), allows for concluding that there are at least four areas where oscilla¬ 
tions frequencies are different. On the other hand, the corresponding maxima 
have not such a regular shape as in the the idealized test example consid¬ 
ered above. For this reason, we introduce the procedure, which extracts the 
irregular regions as 


m = 




Im 


C (a,b)^Rda 

ob 
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where 


C(a, b) = 0 (| w(a, b)\ — L ■ max(|ru(a, 6)|)) 
is the mask with the cut-off threshold L £ [0, 1]. Here 


0(0 = 


i, «e > o 

0, elsewhere. 


is the Heaviside function. 

Thus, the reconstructed features are the results of the inverse continuous 
wavelet transform applied to the regions bounded to be the contours defined 
by the equality | w(a, 6)| — L ■ max(|w(a, 6)|. 

This reconstruction within contours allows for discussing the types of 
dynamics, which are characterized by different frequencies (see Figs l4l3l lower 
panels). In the extracellular signal, as shown in Fig [3j there are oscillations 
with two main frequencies: approximately 6 Hz (within the region marked 
by solid and dashed lines, red, blue and green in the colour online version) 
and 30 Hz (within the black dash-dotted contour). It should be noted, that 
the “fast signal” (thin black line in the lower panel in Figs [3]) is not so strong 
and arises spontaneously at the time about 250 reaching the maximum of 
amplitude in the range of time from 450 to 650. On the other hand, the 
“fast dynamics” within the cell signal (thin line (green in the colour online 
version), FigU lower panel) starts earlier at time about 220 and has a higher 
amplitude. It breaks abruptly and renews at the time f=940. 

Note that our method, which explicitly extracts the main frequency com¬ 
ponents within their regions of switch on and off allows for avoiding their 
mixture with the subthreshold oscillations (so called spikelets), which results 
in the averaged frequency of “fast” signal determined as 118 ± 10 Hz in [lij ]. 

In comparison with the extracellular signal, the “slow” dynamics corre¬ 
sponds to the oscillations of 2 Hz (black dash-dotted line in Fig |4j lower 
panel) with a very small amplitude, which do not have breaks within the 
whole time range of time, to the oscillations of 6 Hz (thick solid (red and 
blue in the colour online version) lines in FigU lower panel)). These oscil¬ 
lations exist within the time intervals from 250 to 400 and from 940 to 1024 
only. 


4. Conclusion and outlook 

In this paper, we have proposed simple computational procedure for the 
inverse continuous wavelet transform that allows for processing of oscillat- 
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Figure 3: The extracellular signal measured in the neocortex layer, its CWT and the partial 
reconstructions of signal’s oscillating features. The plot (upper panel) shows the signal, 
where VFO occurs (the frequency increases in the range of time from 250 to 650). Middle 
panel: Modulus of signal’s CWT, where rectangles mark the regions that contain different 
oscillating components; the contours bound the exact regions used for the reconstruction. 
Lower panel: partial reconstructions of the oscillating components contained in the signal. 
The total time duration corresponds to 4636 ms in the experiment, the maximal amplitude 
is about 0.14 mV 
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Figure 4: The intracellular signal (intrinsic bursting cell (Toj), its CWT and the partial 
reconstructions of intracellular signal’s oscillating features. The plot (upper panel) shows 
the signal, where VFO occurs (the frequency increases in the range of time from 220 
to 650). Middle panel: Modulus of signal’s CWT, where rectangles mark the regions 
that contain different oscillating components; the contours bound the exact regions used 
for the reconstruction. Lower panel: partial reconstructions of the oscillating components 
contained in the signal. The total time duration corresponds to 4636 ms in the experiment, 
the maximal amplitude is about 90 mV 
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ing signals, e.g. the extraction of main frequency components. The modern 
methods of physical research, especially connected with electrophysiology, 
provides data of high complexity, and, correspondingly, their processing re¬ 
quires a development of new approached deeply based of methods of func¬ 
tional analysis (2oJ. The wavelet methods play an important role among 
them. At the same time, they need to be adapted to needs and perceptions 
of computational community. 

Thus, we propose new derivation of the reconstruction formula, which, in 
contrast to the theorem in 16] provides the direct way to a computational 
realization and applications. 

As an example, we have considered the voltage signals measured in neu¬ 
ronal system. A neuronal activity is characterized by the both non-linear 
periodicity and complexity in time and space jT3, EsI, 21, 22 [. Firing pat¬ 


terns could be subdivided into the different frequency bands 21, [19!, 22] and 


connected with the spatial periodicity 0 that provides an important infor¬ 
mation on different types of neurons at pathological processes and a normal 
state. 

We have considered the application of this approach to the real exper¬ 
imental data 19], which describe the spontaneous emergence of very fast 


oscillations at a seizure. The approach considered in the present work allows 
for extracting main periodical components from two signals (extracellular 
and intracellular). The time moments of “switching on” of fast oscillations 
is obtained as well. This information could be useful for the understanding 
of a network structure, in particular, a number of oscillating elements and 
the kind of interconnection between them. 

However, the results are not restricted by the problems of electrophysical 
physiology since the proposed approach and its mathematical and compu¬ 
tational background provide an opportunity to be applied to a much more 
wider class of complex oscillating patterns. 
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Appendix 

MATLAB code, which generates the example represented in Figures [TJ [21 

y„ Determination of the function 
N=512; 

t=linspace(0,1,N); 
fl=cos(20*pi*t); 
f2=exp(-4*t); 

f3=sin(40*pi*t).*ramp(t,t(round(N/3)),t(round(2*N/3))); 
f=f1.*f2+f3; 

°/„ Boundary continuation (extension) 

[te,fe]=fcontin(t,f); 

%Tranform to analytical function 
fa=sAnalytic(fe); 

The corresponding functions, which realize extended equispaced samples 
via the boundary continuation 
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function [tp,fp]=fcontin(t,f); 


n=length(t);tv(:,l)=t;fv(:,l)=f; 

“/.Interpolation into power-two sequences: 

N=2~(ceil(log2(n))); 

ti(:,1)=1inspace(tv(1),tv(end),N); 

fi=interplq(tv,fv,ti); 

“/Boundary continuation 
N=2~(ceil(log2(n))); 
dt=ti(2)-ti(l); 
tp=zeros(2*N,1); 

tp(l:N/2)=linspace(ti(l)-N*dt/2,-dt,N/2); 
tp(N/2+1:1.5*N)=ti; 

tp(l.5*N+1:2*N)=1inspace(ti(end)+dt,ti(end)+N*dt/2,N/2); 
fi=interplq(tv,fv,ti); 

“/„ Boundary reflection 

fp(l:N/2)=flipud(conj(fi(2:N/2+1))); 

fp(N/2+1:1.5*N)=fi; 

fp(1.5*N+l:2*N)=flipud(conj(fi(N/2:N-l))); 

and by the cut off of negative frequencies that forms an analytic function: 
function fa=sAnalytic(f); 

N=length(f); 

“/Cut-off of negative frequencies 
F=fft (f); 

F=[2*F(l:N/2),zeros(l,N/2)] ; 

“/Output: the analytic signal 
fa=ifft(F); 

MATLAB function for the CWT with the Morlet wavelet in the amplitude 
norm, which should be applied to the formed extended sample 

function w=fftMorlet(t,fp,a,omegaO); 

N=length(t); 

“/Fourier transform 
F=fft(fp); 
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nrm=2*pi/(t(end)-t(l)); 

omega_=([(0:(N/2)) (((—N/2)+1):-l)])*nrm; 

7„Convolution 
if a(l)==0 

w(l,:)=fp*exp(-omega0~2/2); 
kl=2; 

else 

kl=l; 

end 

for k=kl:length(a); 

omega_s=a(k)*omega_; 
window=exp(-(omega_s-omegaO)."2/2); 
cnv(k,:)=window.*F; 
w(k,:)=ifft(cnv(k,:)); 

end 

Fig. CD) lower panel) presents abs(w) is obtained as 

omega0=2*pi; 
a=linspace(0,0.2,51) ; 
we=fftMorlet(te,fa,a,omegaO); 
ti=te(0.5*N+1:1.5*N); 
w=we( : , 0 .5*N+1:1.5*N) ; 

The reconstruction procedure is realized in MATLAB as the function 
function iw=invMorlet(t,a,w); 

N=length(t); 
d2t=t(3)-t(l); 

dw( : ,2:N-l) = (w(:,3:end)-w(:,1:end-2))/d2t; 
dw(:,l)=(w(:,2)-w(:,1))/(t(2)—t(1)); 
dw(:,N)=(w(:,N)-w(:,N-1))/(t(N)-t(N-l)); 
iw=imag(trapz(a,dw))/sqrt(2*pi); 
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